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We present a detailed description of the extraction of the highly excited isovector meson spectrum 
on dynamical anisotropic lattices using a new quark-field construction algorithm and a large varia- 
tional basis of operators. With careful operator construction, the combination of these techniques is 
used to identify the continuum spin of extracted states reliably, overcoming the reduced rotational 
symmetry of the cubic lattice. Excited states, states with exotic quantum numbers (0 +_ , 1~ + and 
2" 1 ) and states of high spin are resolved, including, for the first time in a lattice QCD calculation, 
spin-four states. The determinations of the spectrum of isovector mesons and kaons are performed 
on dynamical lattices with two volumes and with pion masses down to ~ 400 MeV, with statistical 
precision typically at or below 1% even for highly excited states. 
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I. INTRODUCTION 

Computing the bound states of QCD is vital if we 
are to claim a complete description of the strong inter- 
actions. Confronting high-precision calculations of the 
spectrum with future experimental measurements will 
test the theoretical framework for such a description rig- 
orously. There has been a resurgence of interest in the 
experimental investigations of the spectrum, notably in 
the charmonium sector where a wealth of high-quality 
data from the _B-factories has challenged our understand- 
ing of spectroscopy. A comprehensive investigation of the 
spectrum of mesons composed of light quarks is the goal 
of the GlueX collaboration, after the 12 GeV upgrade 
of Jefferson Laboratory. Here the aim is to photopro- 
duce mesons, and in particular those with exotic quan- 
tum numbers, as a means of revealing the role of gluonic 
degrees of freedom in the spectrum. 

Lattice calculations offer a method of performing a 
first-principles computation of the spectrum of QCD and 
the calculation of the masses of the lowest-lying states 
has been an important benchmark of lattice studies since 
their inception. However, recently there has been con- 
siderable progress aimed at extracting the spectrum of 
excited states, both for mesons and for baryons. This 
has been accomplished through the use of the varia- 
tional method, employing a large basis of interpolating 
operators satisfying the symmetries allowed by the cubic 
lattice[THl]. In a series of recent papers, we have applied 
this methodology to the extraction of the meson spec- 
trum [5], and the radiative transitions between excited 
and low-lying meson states [6J [7] . The first studies were 
performed in the quenched approximation to QCD, for 
mesons composed of the heavier charm quark and its anti- 
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quark, a system which is computationally less demanding 
yet for which there is a wealth of high-quality experimen- 
tal data. In this paper, we investigate the spectrum of 
mesons for quark masses below the strange quark mass, 
going down to pion masses of around 400 MeV, and ex- 
pand on our earlier letter [5] focused on the spectrum of 
mesons in full QCD with three mass-degenerate quark 
flavors. 

Several previous studies of the spectrum have focussed 
on obtaining precision results for ground state masses 
through controlling systematic errors [HHHJ- Others have 
aimed at extracting the masses of some of the excited 
states [T5HTU] . The extraction of excited state masses 
is more difficult owing to the decrease in the signal-to- 
noise ratio with increasing time as we move progressively 
higher in the spectrum. To circumvent this difficulty, we 
use anisotropic lattices, with finer temporal than spatial 
lattice spacing, enabling the behavior of the Euclidean- 
space correlation functions to be examined at small tem- 
poral separations. 

The (hyper-) cubic lattice does not possess the full ro- 
tational symmetry of the continuum. Thus in a lattice 
calculation, states at rest are classified not according to 
the spin (J,J Z ), but rather according to the irreducible 
representations (irreps) of a cube; for states of higher 
spin, the different continuum degrees of freedom are dis- 
tributed across several lattice irreps. In this study, we use 
a large basis of interpolating operators, decomposed into 
their lattice irreps, that enable us to explore all J PC up to 
spin 4, except for the exotic 4" 1 , with as many as 26 op- 
erators in a given symmetry channel. Our ability to cal- 
culate correlation functions efficiently for such a large ba- 
sis of operators relies on a new method, "distillation" [20] , 
for the construction of quark-antiquark operators, includ- 
ing those with non-local construction. 

As the lattice spacing approaches zero, full rotational 
symmetry is restored and thus in principle the spins can 
be identified through the emergence of energy degenera- 
cies between different irreps. The increasing density of 
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states in each irrep as we rise in the spectrum makes the 
identification of such degeneracies a challenging task, and 
thus the assignment of the continuum spins has been a 
formidable barrier to the interpretation of lattice calcula- 
tions. In this paper, we adopt a very different approach, 
in which by judiciously constructing operators so as to 
have a known continuum behavior, the spins of the ex- 
cited states are determined and thus the barrier imposed 
by the reduced cubic symmetry of the lattice is overcome. 

A particularly interesting class of mesons are those 
with exotic J PC , such as the 1 h channel. Extracting 
clean signals for even the lightest state with these quan- 
tum numbers has proven difficult [HI I2TH26] , with statis- 
tical noise levels typically being significantly higher than 
for other states. There is also a requirement to determine 
that any such extracted state is, in fact, the exotic spin 
1 and not a non-exotic 4 h state which would live in the 
same lattice irreducible representation. We found that 
we can extract information about exotic state masses at 
the same level of precision as excited non-exotic states 
(< 1%), with the spins clearly identified. 

First results, exploiting the full planopy of anisotropic 
lattices, "distillation" for efficient computation of the in- 
terpolating operators, and the identification of the con- 
tinuum spins, have been presented for the case of three 
degenerate "strange" quarks [5] ; the spectrum of excited 
states, including those with high spin, was extracted with 
confidence. 

In this paper, we expand on the earlier work (at 
ss 700 MeV) to include calculations in full QCD 
both with three degenerate flavors of quark, and with 
a strange quark and two light quarks (Nf — 2 1), cor- 
responding to pion masses ~ (520,440,400) MeV. 
Furthermore, we perform calculations at two spatial vol- 
umes, enabling us to seek possible finite-volume effects, 
and potentially the presence of multi-hadron states in the 
calculated spectrum. Precise extraction of such states is 
required in order to carry out analyses of resonances in 
the manner suggested by Luscher [27] . 

The structure of the paper is as follows. We begin in 
Section [TT] by presenting details of the QCD gauge-field 
configurations to be used. We describe the technology of 
two-point correlator measurement using distillation on 
dynamical gauge-field configurations in Section |III| In 
Section |IV| the construction of a set of composite QCD 
operators suitable for use in the extraction of the me- 
son spectrum is outlined, along with a description of the 
procedure used to make these operators transform irre- 
ducibly under the limited rotations allowed on a cubic lat- 
tice. In Section|V]we present the methodology utilised to 
extract meson spectral information from correlation func- 
tions. Section |VI| discusses the possibility of using the 
information embedded in vacuum-opcrator-state matrix 
elements to determine the spin of a meson, overcoming 
ambiguities introduced by the reduced rotational sym- 
metry of a cubic lattice. The stability of the extracted 
spectral quantities to changes in the details of the cor- 
relator analysis is considered in Section |VII| where it is 
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TABLE I: The lattice data sets and propagators used in this 
paper. The lattice size and number of configurations are 
listed, as well as the number of time-sources and the number 
of distillation vectors N vccs . The total number of inversions 
for each quark mass = 4 x AT tijros x N vocs x -/V c f gs is shown. 



seen that we can extract results of considerable robust- 
ness. In Section VIII we report meson spectrum results 



extracted from calculations at four pion masses and two 
different lattice volumes. Results for isovector mesons, 
kaons and the connected part of ss ( "strangeonium" ) 
are shown. Section |IX| considers the apparent absence 
of multi-particle states within our extracted spectra. Fi- 
nally, in Section [X] we summarise our observations and 
suggest future directions. 



II. LATTICE GAUGE FIELDS 

In Euclidean space, excited state correlation functions 
decay faster than the ground state, and at large times 
are swamped by the signals of lower states, thus com- 
plicating the resolution of excited states. To ameliorate 
this problem we have adopted a dynamical anisotropic 
lattice formulation whereby the temporal extent is dis- 
cretized with a finer lattice spacing than in the spatial 
directions [35J HH] ; this has proven crucial to obtain the 
results shown in this paper. This method avoids the 
computational cost that would come from reducing the 
spacing in all directions. Improved gauge and fermion 
actions are used, corresponding to two light dynamical 
quarks and one strange dynamical quark. Details de- 
scribing the formulation of the actions as well as the 
techniques used to determine the anisotropy parameters 
can be found in Refs. [25J EH] ■ The lattices have a spatial 
lattice spacing a s ~ 0.12 fm with a temporal lattice spac- 
ing 3.5 times smaller corresponding to a temporal scale 
a" 1 - 5.6 GeV. 

Previous work [5J showed results using the three-flavor 
degenerate quark-mass dataset corresponding to bare 
light and strange quark masses atini = atm s = —0.0743 
and lattice size 16 3 x 128. The pion mass (degenerate 
with the kaon and rj masses) is roughly 700 MeV. In 
this work, results are extended to lighter masses a t mi = 
(-0.0808, -0.0830, -0.0840) and a t m s = -0.0743 corre- 
sponding to 2 + 1 flavors of dynamical quarks, and lattice 
sizes 16 3 x 128 as well as 20 3 x 128. These fully dynam- 
ical datasets, described in more detail in Table [IJ allow 
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for some investigations of the quark mass dependence as 
well as finite volume dependence of the spectrum. 

The lattice scale, quoted above, is determined in the 
physical quark mass limit using the f2 baryon mass (de- 
noted by atmn). As noted, these calculations are carried 
out away from this limit. To facilitate comparisons of the 
spectrum at different quark masses, the ratio of hadron 
masses with the f2 baryon mass is used to remove the 
explicit scale dependence [2§] . 



III. CORRELATOR CONSTRUCTION 

The determination of the excited spectrum proceeds 
from the calculation of correlation functions between a 
basis of hermitian creation and annihilation operators O 
at Euclidean times and t, 



ci i (t) = (o|o,(t)o J (o)|o). 



Inserting a complete set of eigenstates of the Hamilto- 
nian, such that H\n) = E n \n), this correlation function 
decomposes into a sum of contributions from all states 
in the spectrum with the same quantum numbers as the 
source operators, 
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C ii(*)=E^ n < l^l n ^ n l°il°) e ~ 



E„t 



where the discrete character of the spectrum follows be- 
cause the calculation is performed in finite-volume. In 
order to measure energies of low-lying states, it is crucial 
to construct operators that overlap predominantly with 
light modes. 

Smearing is a well-established means to improve op- 
erator overlap, whereby a smoothing function is applied 
to the quark fields used in the creation operators. This 
smoothing function should effectively remove noisy short- 
range modes which should not make a significant contri- 
bution to the low-energy correlation functions. The Ja- 
cobi method [3D] uses the gauge-covariant second-order 
three-dimensional lattice Laplacian operator 
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where the gauge fields, U may be constructed from an ap- 
propriate covariant gauge-field-smearing algorithm [31 j . 
To suppress high-energy modes of V 2 , this operator is 
exponentiated, exp (cV 2 ), with some smearing weight a. 
The resulting smoothed operator is then applied to the 
quark fields ip. 

The suppression of the high energy modes of the Jacobi 
smearing operator exp(crV 2 ) means that only a small 
number of modes contribute significantly to the construc- 
tion of the smeared quark fields, ip. As suggested in 
Ref. |20) . this smearing function can be replaced with 



a low-rank approximation. The "distillation" operator 
defines a smearing function 



A? 



= E^ fc) (*)4 fc)t (*). 



(i) 



k=l 



where £i are a finite number, N, of eigenvectors of V 2 
evaluated on the background of the spatial gauge-fields 
of time-slice t, once the eigenvectors have been sorted 
by eigenvalue. This is the projection operator into the 
subspace spanned by these eigenmodes, so D 2 = □. 
This smearing function is used in the construction of 

isovector meson operators of the form ipTip, where T acts 
in spin and color as well as coordinate space. Applying 
the distillation operator □ onto each quark field, ip = 
dip, the creation operators at zero three-momentum arc 
written as 

Oi(t) = Mt)n xy {t) ■ r yz {t) ■ u zw (t)ip w (t), 

where there is an implied volume summation over re- 
peated spatial indices. In a shorthand notation the cor- 
relation function can be written as 



After evaluating the quark-field path-integral and insert- 
ing the outer-product definition of the distillation oper- 
ator □ from Eq. [T] the correlator can be written 



dj{t) = Tr $ J (t)r(t,0)$ t (0)r(0,t) 

where 

^ /3 (o) = C t (o) [r(o)] Q ^(o), 

encodes the structure of the operator and 

7-afl(*,0)=^(i)M- /s 1 (* > 0)C(0), 



is the "permabulator" , with M the lattice representation 
of the Dirac operator and where the quark spin indices, 
a, j3 of $ and r have been explicitly written. 4> has a well- 
defined momentum, while there is no explicit momentum 
projection in the definition of r. The $ and r are square 
matrices of dimension NN& where N a = 4 are the num- 
ber of spin components in a lattice Dirac spinor. Con- 
struction of the t require NN a inversions of the fermion 
matrix to compute all elements. These matrices are small 
compared to the dimension of the space of quark fields. 
Once the r have been computed and stored, the correla- 
tion of any source and sink operators can be computed 
a posteriori. The method straightforwardly extends to 
the determination of multi-hadron two-point correlation 
functions [35] as well as three-point functions |20) . 
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IV. CONSTRUCTION OF MESON OPERATORS 
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Meson spectral information will follow from analysis of 
two-point correlators featuring a large basis of composite 
QCD operators having mesonic quantum numbers. The 
simplest such operators are color-singlet local fermion bi- 
linears, il>i a (x,t)Y a pipip(x,t) where the quantum num- 
bers are determined by the choice of gamma matrix, T. 
In distillation [2U] the quark fields, tp are replaced by 
the smeared quark fields ip but the rotationally symmet- 
ric nature of the smearing does not change the quan- 
tum numbers of the bilinear operators. These simple 
local operators are extremely limited in that they allow 
access only to the set J PC = 0~+, 0++, 1 — , 1++, 1+" 
and they do not offer significant redundancy within any 
J PC . In order to consider higher spins, exotic J PC and 
to produce multiple operators within a given symmetry 
channel, one must consider extending to the use of non- 
local operators [I3HT8] . Our approach is to use spatially- 
directed gauge-covariant derivatives within a fermion bi- 
linear, that is to construct operators of essential structure 



.ip(x*,t) 



where D 



and where spin and color indices 
are suppressed for clarity. The use of the "forward- 
backward" derivative, D is not strictly necessary at zero 
momentum (projected by the sum over spatial sites), the 
only case we consider here, but it does somewhat simplify 
the construction of eigen-operators of charge-conjugation 
as will be discussed below. 

With the continuum SO (3) rotational symmetry it is 
straightforward to produce operators of this type that 
are of definite spin, parity and charge-conjugation at zero 
momentum. This follows from forming a circular basis of 
the cartesian-vector-like derivatives and gamma matri- 
ces, Di,ji, 757j, 7o7i, tijk'Yjlk-, e.g. 



u m =- 1 - 75 

D m=0 = 1 D 

D m =+l = ~ 



V2 



i — > 



■iD, 



TABLE II: Gamma matrix naming scheme. 



At the two-derivative level we adopt the convention of 
first coupling the two derivatives to a definite spin, Jd, 
then coupling with the vector-like gamma matrix (if any) 
as 



[2] xJ,M 



E 

m 3 ,m D 



(l,m 3 ; J D ,m D \j,M) 



x (l,mi; l,m 2 \JD,m D ) ipY rrl3 D mi D m2 ip. 

It is worth noting here that while the Clebsch-Gordan for 
1 ® 1 — ¥ 1 is antisymmetric and D TOl D m2 appears to 
be symmetric, there are non-zero operators with Jd = 1 
because the gauge-covariant derivatives do not commute 
with each other in QCD. Rather these "commutator" 
operators are proportional to the gluonic field-strength 
tensor which does not vanish on non-trivial gluonic field 
configurations. 

At the three-derivative level we need to choose a con- 
vention for the ordering in which we couple the deriva- 
tives. A natural choice comes from insisting the op- 
erators have definite charge-conjugation symmetry. By 
exchanging the quark and anti-quark fields, the charge- 
conjugation operation effectively acts as a transpose of 
the operators between the quark fields - for three deriva- 
tives then, one ensures definite charge conjugation by 
coupling the outermost derivatives together first since 
this gives them a definite exchange symmetry (even for 
J13 = 0,2, odd for J 13 = 1). This simple formulation is 
possible because we have used the "forward-backward" 
derivatives, D . 



^ (l, m 4 ; J D ,m D \j 1 M)(l, m 2 ; J13, TOj. 3 | J D , m D ) 

mi ,771-2 , 
7713 ,7774 , 
77113, 771 

x (l,mi; l,m 3 |Ji3, mia)ipT mi D mi D „ l2 D rri3 ip. 



Once expressed in this basis, which transforms like spin- 
1, operators of definite spin can be constructed using the 
standard SO(3) Clebsch-Gordan coefficients. For exam- 
ple, with a vector-like gamma matrix and one covariant 
derivative, operators of J = 0, 1, 2 can be formed 



mi ,m2 



(l,mi;l,m 2 |J,M) ipT mi D m2 ^. 



Clearly this procedure can be extended to as many 
covariant derivatives as one wishes. In this paper we 
will use operators with up to three derivatives providing 
access to all J PC with J < 4. 1 

The operators as formed are eigenstates of charge- 
conjugation in the case that the ip and ip fields are of 
the same flavour. In the case that the fields are degen- 
erate but not identical (e.g. the u and d quarks in our 



The choice of Y plays a role in setting the parity and 
charge-conjugation quantum numbers of the operator - 
our naming scheme for these matrices is given in Table 



1 except the exotic 4+ which requires a minimum of four deriva- 
tives. 
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TABLE IV: Number of operators in each lattice irrep A 
TABLE III: Continuum spins subduced into lattice irreps using all operators with up to three derivatives. 
A(dim). 



calculation), the C-parity is trivially generalised to G- 
parity. For kaons, where the light and strange quarks are 
not degenerate, there is no C-parity or any generalisation 
of it. In this case the symmetry channels are labelled by 
J p and operators of both C can be used. 



A. Subduction into lattice irreps. 

In lattice QCD calculations the theory is discretized on 
a four-dimensional hypercubic Euclidean grid. The full 
three-dimensional rotational symmetry that classifies en- 
ergy eigenstates in the continuum is hence reduced to the 
symmetry group of a cube (the cubic symmetry group, or 
equivalently the octahedral group) . Instead of the infinite 
number of irreducible representations labelled by spin </, 
the single-cover cubic group relevant for integer spin has 
only five irreducible representations {irreps): Ax, T%, T%, 
E, A 2 . The distribution of the various M components of 
a spin- J meson into the lattice irreps is known as sub- 
duction, the result of which is shown in Table III 



To be of any use in lattice computations, the contin- 
uum operators described above must be subduced into 
lattice irreps. Noting that each class of operator is closed 
under rotations, the subductions can be performed using 
known linear combinations of the M components for each 
J: 



O 



J2 s if (rx^J) J > M = ^s 



J,MqJ,M 



A, A 



A I 



A I 



where A is the "row" of the irrep (1 . . . dim(A)). Note 
that, although oj/i can have an overlap with all spins 



contained within A (as listed in Table III for J < 4) 



it still carries the memory of the J from which it was 



subduced, a feature we exploit in Section [VI] The sub- 
duction coefficients, S ,JW 

^J,M C J,M* 



End,JW oJ 
AI °A,A °A',\> 



. > , form an orthogonal matrix, 

<5a,A'<5a,A' , and this fixes their nor- 
malisation. 

The subduction coefficients can be constructed in a 
number of different ways and here we give a simple 
derivation. More details and an alternative method us- 
ing a group-theoretic projection formula can be found in 
Appendix [A] 

The simplest case is the subduction of the J = oper- 
ator; from Table III this only subduces into the A\ irrep 



and so trivially we have 1 = 1. The J — 1 opera- 
tor is also relatively straightforward, only subducing into 
the the T\ irrep with subduction coefficients shown in 



M 
Ti,A 



= <5a,2-a/, where the shift 
dim(A). 



Appendix A Note that S\ 
by 2 placesA in the range 1 . 

Subduction coefficients for all higher spins can be con- 
structed by iteration, starting from the J — and J = 1 
coefficients and using 



S 



J,M 
A, A 



= N 



oJi,Mi C J 2 ,M 2 



A 1: A 2 Mi,M ; 

A Ai A 



Ai,A! °A 2 ,A 2 



°( A Ai Xo) {Jl ' Mi;J2 > M2lJ ' M) 



Here (Ji, M x \ J 2 , M 2 \J, M) is the usual SO(3) Clebsch- 

Gordan coefficient for J x (g) J 2 — > J and c( ^ ^ 2 \ 

^ A Ai A 2 ' 

is the octahedral group Clebsch-Gordan coefficient for 
Ai (£> A2 — > A. N is a normalisation factor, fixed by 
the requirement that the subduction coefficients form an 
orthogonal matrix as discussed above. We give explicit 
values for the subduction coefficients up to J = 4 in 
Appendix [X] 

In Table |IV| we show the number of operators we have 
in each lattice irrep, i.e. using all operators with up to 
three derivatives. We have performed extensive tests of 
this operator set to check that two-point correlators hav- 
ing operators in differing irreps at source and sink are 
consistent with zero and that similarly within an irrep, 
correlators of differing rows at source and sink are consis- 
tent with zero. All such "orthogonality tests" are passed 
in explicit calculation. 

For our final spectral extractions, we form a correlator 
matrix in a given irrep A and average over equivalent 
rows, A, 



dim(A) 

Cij = dim(A) X] = dimfA^ ^^1^(A)A^3(A)a1^)j 

A=l A 



where i,j labels the different operator constructions 
within the irrep A. 
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V. CORRELATOR ANALYSIS 

The variational method for spectral extraction [33l 134] . 
which takes advantage of a redundancy of operators 
within a given symmetry channel, is now in common 
usage [SJ U31 [T5HT7] . This method finds the best (in a 
variational sense) linear combination of operators within 
a finite basis for each state in the spectrum. Mathemat- 
ically it boils down to the solution of a linear system of 
generalised eigenvalue type: 

c{ty\t) = x n (t)c(t )v n (t) (2) 

where A tt (to) = 1 and where there is an orthogonality 
condition on the eigenvectors of different states (n, n'), 
v n — ^n,n'- As discussed in Ref. [3] this or- 
thogonality condition is very powerful in extracting near 
degenerate states which would be difficult to distinguish 
by their time dependence alone. 

In our particular implementation of this method, equa- 
tion [5] is solved for eigenvalues A n and eigenvectors v n , 
independently on each timeslice, t. Ensuring the same 
ordering of states between timeslices requires some care 
owing to the high degree of mass degeneracy in the me- 
son spectrum. Rather than the obvious ordering by size 
of eigenvalue which might fluctuate timeslice-by-timcslicc 
for nearby masses, we associate states between timeslices 
using the similarity of their eigenvectors. We choose a 
reference timeslice on which reference eigenvectors are de- 
fined, i/" ef = v n (t Te {), and compare eigenvectors on other 

timeslices by finding the maximum value of v^J C(to)v n 
which associates a state n with a reference state n'. Us- 
ing this procedure we observe essentially no "flipping" 
between states in cither the principal correlators, A n (i) 
or the eigenvectors v n (t), as functions of t. 

Within finite-volume field theory, which has only dis- 
crete eigenstates, any two-point correlator can be ex- 
pressed as a spectral decomposition, 

where this is an approximation valid providing f < I ( , 
the temporal length of the box. The "overlap factors", 
Z" = (n|Oi|0) are related to the eigenvectors by Z" = 
V2^"e m " 4 °/ 2 v^Cjiito). The state masses follow from 
fitting the principal correlators, A n (i), which for large 
times should tend to e~ mn ^~ to ^. In practice we allow a 
second exponential in the fit form and use even relatively 
low timeslices in order to stabilise the fit. The fit function 
is 

A n (i) = (1 - A n )e-" l " (t - to) + A n e- m '« {t - to \ (4) 

where the fit parameters are m n , m' a and A n . Typical fits 
for a set of excited states within an irrep are shown in 
figure [l]whe re we plot the principal correlator with the 
dominant time-dependence due to state n divided out. In 
such a plot one would see a horizontal line of value 1.0 in 
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FIG. 1: Principal correlator fits according to eqnQ. 
Eight states from the Tf" irrep (743, 16 3 ). Plotted are 
X n (t) ■ e m " ( '" to) data and the fit for to = 8. Data used 
in the fit are shown in black, while points excluded from the 
fit are in grey. 



the case that a single exponential dominates the fit and 
clearly the data shows fiat behaviour for t > to. 

Empirically we find that the size of the second expo- 
nential term decreases rapidly as one increases to- Fur- 
ther we find, in agreement with the perturbative analysis 
of Ref. [35 , that for large to values the m' a extracted are 
larger than the value of m n=( ji m (c) , the largest "first" ex- 
ponential mass extracted. At smaller to values this is not 
necessarily true and is indicative of forcing an incorrect 
orthogonality as discussed below. The values of A n and 
m' n are not used elsewhere in the analysis. 

From the spectral decomposition of the correlator, 
equation [3j it is clear that there should in fact be no time 
dependence in the eigenvectors, while our independent 
solution of the generalised eigenvalue problem as a func- 
tion of t has allowed there to be. The time-independent 
overlap factors, iT", which will be used later to identify 
the spin of extracted states, follow from fitting the Z™(t), 
obtained from the eigenvectors, with a constant or a con- 
stant plus an exponential (in the spirit of the perturbative 
corrections outlined in |35|). 

The importance of choosing an appropriately large 
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FIG. 2: Normalised corrrelation matrix (C\j / \jCuCjj) on 
timeslice 5 in the T x irrep (743, 16 3 ). Operators are ordered 
such that those subduced from spin 1 appear first followed by 
spin 3 then spin 4. 

value of t was emphasised in Ref. [5] . In this paper we 
will follow the "reconstruction" scheme outlined therein 
in the selection of to- In short, the masses, m a , extracted 
from fits to the principal correlators and the Z" extracted 
from the eigenvectors on a single timeslice are used in 
equation [3] to "reconstruct" the correlator matrix. This 
reconstructed matrix is compared to the original data 
for all t > to with the degree of agreement indicating the 
acceptability of the spectral description. The description 
generally improves as one increases to until at some point 
the increase in statistical noise prevents further improve- 
ment. In particular see figure 6 in Ref. [5. where the 
effect of choosing t too small is clearly seen. Forcing 
the dim(C)-state orthogonality, v ml< C(to) v — o n m , m 
a situation where accurate description of C(to) requires 
more than dim(C) states leads to a poor description of 
the correlator matrix at times t > to- The reconstruction 
procedure gives a guide to the minimal to for which the 
correlator matrix is well described by the variational so- 
lution. The sensitivity of extracted spectral quantities to 
the value of to used will be discussed in detail in section 
VII A| but in short it is usually necessary for us to use 
to >7. 

The reconfit2 code used for variational analysis is 
available within the adat suite [55] . 



VI. DETERMINING THE SPIN OF A STATE 

In principle the most rigourous method to determine 
the spin of a state is to perform the extraction of the 
spectrum for each lattice irrep at successively finer lat- 
tice spacings, and then to extrapolate the energies in 
each irrep to the continuum limit. There one expects 
to see degeneracies emerge according to the pattern of 
subduction, free of splittings arising from the discreti- 
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0.612(1) 1.113(7) 1.188(4) 1.207(5) 1.356(10) 1.560(7) 1.584(12) 1.620(12)1.648(23) 

L1LJJ 

1.201(4) 1.204(4) 1.554(25) 1.576(16) 1.626(8) 

IJj 

1.190(4) 1.565(25) 1.577(21) 1.603(26) 

FIG. 3: Overlaps, Z, of a selection of operators onto states 
labelled by m/mn in each lattice irrep, A {743, 16 3 ). Z's 
are normalised so that the largest value across all states is 
equal to 1. Lighter area at the head of each bar represents 
the one sigma statistical uncertainly. 



sation effects. Thus, for example, a spin-3 state would 
appear as degenerate energies within the A2 , Xi and T2 
irreps. This procedure has been successfully applied to 
identify a number of low-lying states in the calculation 
of the glueball spectrum within pure SU(3) Yang-Mills 
theory [37]. 

There are two reasons why this technique is not cur- 
rently practical for the QCD meson spectrum. Firstly, 
the procedure relies on a series of calculations on pro- 
gressively finer lattices, and hence at increasing computa- 
tional cost. Secondly, the continuum spectrum, classified 
according to the continuum quantum numbers, exhibits 
a high degree of degeneracy; when classified according 
to the symmetries of the cube, the degree of degener- 
acy is vastly magnified. Identification of degeneracies 
between irreps would require a statistical precision far 
beyond even that of the high-quality data presented here, 
as seen in Figure [TO] and subsequent figures. 

To alleviate these difficulties it would be useful to have 
a spin-identification method that is effective when us- 
ing data obtained at only a single lattice spacing. Ob- 
viously this lattice spacing should be fine enough that 
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FIG. 4: A selection of Z values across irreps A for states suspected of being J = 2, 3, 4 (743, 16 3 ). Left to right the operators 
are (ai x D^ =i ) J=2 ,{p x £>£L) J = 2 , (p 2 x D [ f =2 ) J =\ (a x D l ^ =2>J=2 ) J = 2 , (6 x D^, =1 , J=a ) J = 2 , (ai X 7?^^) J = 2 , (p x 
Df =2 ) J = 3 , (P2 x D?i 2 ) J = 3 , (a x Df i3=2>J=3 ) J = s , (ai x ^ 3=2 , J=3 ) J = 3 , (oi x ^,=3,^) J=3 > ^ x ^=1,^) J=3 and fa * 



rotation symmetry has been restored to a sufficient de- 
gree in order that it be describing QCD. The mass de- 
generacy complications outlined above suggest that any 
alternative method needs to use state information beyond 
just the mass. Our suggestion is to consider the values 
of the vacuum-to-state matrix elements, or "overlaps" 
((n|O|0)) of our carefully constructed subduced opera- 
tors. 



The operators constructed in Section IV transform ir- 



reducibly under the allowed cubic rotations, that is they 
faithfully respect the symmetries of the lattice. How- 
ever it is also clear from the method of construction that 
each operator carries a "memory" of the contin- 

uum spin, J, from which it was subduced. If our lat- 
tice is reasonably close to restoring rotational symmetry 
we would expect an operator subduced from spin J to 
overlap strongly only onto states of continuum spin J. 
In fact this is clearly apparent even at the level of the 
correlator matrix as seen in figure [2j Here the correlator 
matrix for T x is observed to be approximately block di- 
agonal when the operators are ordered according to the 
spin from which they were subduced. The elements out- 
side the diagonal blocks are smallest when the operators 
feature zero, one or two derivatives and are somewhat 
larger for three-derivative operators, possible reasons for 
this will be discussed later in this section. 

The effect is seen even more strongly at the level of 
individual states, where the "overlaps", Z" — (n|Oi|0) 
for a given state show a clear preference for overlap onto 
only operators of a single spin. In figure [3] we show the 
overlaps for a set of low-lying states in the A irreps of 
the 743 16 3 calculation, the mass spectrum being shown 
in the second pane of figure [lO] While we show only a 
subset of the operators for clarity, the same pattern is 
observed for the full operator set. 

In fact we can be more quantitative in our analy- 
sis and compare the overlaps obtained in different ir- 
reps. In the continuum our operators are of definite spin 
such that (O|0 J < M |J',Af) = Z^5,j tJ ,5 MM , and there- 



fore (0|O^lj J', M) = S^Z^Sjj, so that only the spin 
J states will contribute, and not any of the other spins 
present in the irrep A. Inserting a complete set of me- 
son states between the operators in the correlator and 
using the fact that the subduction coefficients form an 
orthogonal matrix, Y^m^a^^A' X 



S A,\ S A',y = "A, A' "a, A' 1 we ob- 
tain terms in the correlator spectral decomposition pro- 
portional to Z^*Z^; these terms do not depend upon 
which A we have subduced into, up to discretisation un- 
certainties as described below. Hence, for example, a 
J = 3 meson created by a [J = 3] operator will have 
the same Z value in each of the A 2 ,Ti,T 2 irreps. This 
suggests that we compare the independently obtained Z- 
values in each irrep. In figure [4] we show the extracted Z 
values for states suspected of being spin 2,3 and 4 across 
the A irreps. 

As seen in figure [3J at finite lattice spacing there are 
deviations from exact equality. Some discretisation ef- 
fects scale with positive powers of the lattice spacing, 
such as the effect of using finite differences to represent 
derivatives. There are no dimension-five operators made 
of quark bilinears that respect the symmetries of lattice 
actions based on the Wilson formalism and that do not 
also transform trivially under the continuum group of 
spatial rotations. Thus, rotational symmetry breaking 
terms do not appear until 0(a 2 ). This argument holds 
even though the action used in this work describes an 
anisotropic lattice. As a result, we expect the rotational 
breaking between lattice irreps to be suppressed in both 
the spectrum as well as for the wave- function overlaps. 
On the other hand, renormalisation mixing of high mass- 
dimension operators with lower mass dimension opera- 
tors can give rise to effects scaling with negative powers 
of the lattice spacing and these we would expect to cause 
more trouble when Z values are to be used to determine 
spin. In practice we do see the largest discrepancies for 
operators with three derivatives, but even here the ef- 
fects are not so large as to prevent use of the method. 
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FIG. 5: Fit to the four subduced principal correlators of a 4 ++ meson using a common mass (743, 16 3 ). Plotted is A(t)-e 
Grey points not included in the fit. 
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We suggest that it is our use of smoothed fields 2 that has 
rendered these mixings relatively small, sensitive as they 
are to high-energy physics which has been filtered out. 

In summary, we have demonstrated that the Z values 
of carefully constructed subduced operators can be used 
to identify the continuum spin of states extracted in ex- 
plicit computation, as least on the lattices we have used. 

Given that this is possible, suppose we confidently 
identify the components of a spin- J meson subduced 
across multiple irreps; what then should we quote as our 
best estimate of the mass of the state? The mass de- 
termined from fits to principal correlators in each irrep 
can differ slightly due to unavoidable discretisation ef- 
fects and avoidable fitting fluctuations (such as variations 
in fitting time-region). In practice we have found that 
variations due to changes in the fitting of principal cor- 
relators are typically much larger than any discretisation 
differences and we propose a simple scheme to minimise 
these in a final quoted mass. Rather than averaging the 
masses from independent fits to multiple principal corre- 
lators, we perform a joint fit to the principal correlators 
with the mass being common. We allow a differing sec- 
ond exponential in each principal correlator so that the 
fit parameters are m n ,{m n A } and {^n}- These fits are 
typically very successful with correlated xV-^dof close to 
1. An example for the case of 4 ++ components identified 
in A+ + ,T+ + ,T+ + ,E++ is shown in figure |J When we 
present our final, spin-assigned spectra it is the results of 
such fits that we show. 



VII. STABILITY OF SPECTRUM 
EXTRACTION 

In this section we consider to what extent the extracted 
spectrum changes as we vary details of the calculation, 
"keeping the physics constant" . Variations to be consid- 
ered are the specific reference timeslice, to, used in the 
variational analysis, the set of meson operators used and 
the number of distillation vectors. We will use the T x 
irrep in the 743, 16 3 dataset to demonstrate our findings. 



2 All gauge links are stout-smeared and the distilled quark fields 
are effectively low-momentum filtered 



A. Variational analysis and to 

Our fitting methodology was described in Section [V] 
where reconstruction of the correlator was used to guide 
us to an appropriate value of to. As seen in Figure [6j 
f° r ^0 ^ 6, the low-lying mass spectrum is rather stable 
under variations of to. This appears to be mostly due to 
the inclusion of a second exponential term in equation [4] 
which is able to absorb much of the effect of other states 
"leaking" into this principal correlator through use of an 
inaccurate orthogonality. The contribution of this second 
exponential typically falls rapidly with increasing t Q both 
by having a smaller A and a larger ml . 

Overlaps, Z\ = (n|O,|0), can show more of a sensitiv- 
ity to to values being too low as one might expect given 
the argument of an incorrect orthogonality in the gen- 
eralised eigenvector space at small to . In Figure [7] we 
present overlaps of various J = 1 states onto an operator 
subduced from J = 1 and the overlap of an extracted 
J = 4 state onto the only J PC = 4 operator in our 
basis. Clearly in the J = 4 case one only extracts a stable 
Z for large £0 , which is likely due to heavier J = 4 states 
only here becoming negligible contributions to C(tn). 

We note that one may fit the extracted Z{to) with 
either a constant or a constant plus an exponential as 
shown in figure [7] Since the data between to values are 
strongly correlated, the statistical uncertainty is not sig- 
nificantly decreased by this procedure, but it does seem 
to "average out" some of the fluctuations in fitting at 
each to- 
In summary it appears that variational fitting is 
reliable provided to is "large enough". Using two- 
exponential fits in principal correlators we observe rel- 
atively small to dependence of masses, but more signif- 
icant dependence for the Z values which we require for 
spin-identification. 



B. Changing the operator basis 

Our approach, as described in Section |IV| is to con- 
struct a variational correlator matrix featuring all opera- 
tors available to us in a given irrep at up to three deriva- 
tives. Here we consider the effect on the extracted spec- 
trum of reducing the size of this operator basis. In the 
plots that follow we will use the color-coding described 
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FIG. 6: Extracted 7\ mass spectrum as a function of to- 
Horizontal bands to guide the eye. 



in Table |V| to indicate spin-assignments. 

In the first four columns of Figure [SJ the operator basis 
is reduced by discarding operators, first discarding those 
whose diagonal correlators are noisiest and then, in col- 
umn (d), discarding all three-derivative operators except 
the one J = 4 operator. We observe that it is only high 
in the spectrum that any change takes place and that 
there it is only at the level of statistical fluctuations. In 



(a) 




to/a, 



lb) 




(d) the first excited J = 3 state is no longer extracted 
- this is most likely due to the basis only retaining two 
J = 3 operators. 

In the next four columns we deliberately discard oper- 
ators we believe are important. In column (e), all opera- 
tors featuring commutators of derivatives are discarded. 
As discussed later in this manuscript these operators 
may have good overlap onto states containing an excited 
gluonic field. The state at m/rriQ ~ 1.35 is observed, 
using the full basis of operators, to have large overlap 
onto a commutator operator, (n x £)?]_■, ) J=1 along with 
smaller overlaps onto non-commutator operators such as 
p. In column (e), the commutator operators being dis- 
carded, we find that this state is less cleanly extracted 
when it can only be produced through its suppressed non- 
commutator overlaps. Otherwise the spectrum in (e) is 
rather similar to that obtained with the full operator ba- 
sis. 

In column (f) all operators subduced from spin J = 3 
are discarded. The observed spectrum is almost identical 
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FIG. 7: Extracted overlaps as a function of to. Fits to a 
constant or to a constant plus a decaying exponential shown 
by the coloured regions. 
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FIG. 8: Extracted T" mass spectrum for various operator 
bases, (a)-(d) are "reasonable" operator bases, (e)-(h) discard 
important operators. 

(a) Full basis (dim = 26), (b) Full basis less two noisiest O j=1 
and noisiest J=3 (dim = 23), (c) Full less four noisiest C J=1 
and two noisiest C J=3 (dim = 20), (d) No three-derivative 
operators except 0' 1=i (dim = 13), (e) No operators with 
commutators of derivatives (dim = 15), (f) No J ~ i (dim = 
20), (g) No Q J - iA (dim = 19), (h) No J=i (dim = 25). 
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FIG. 9: Extracted 7\ mass spectrum as a function of number of distillation vectors. Green points and inset show the effect 
of using a modified smearing operator a w = J2n=i w (^n)£ n £ n * with w(X„) = e~ a x ™' 4 . 



to that with the full basis with the expected exception 
of the previously identified J = 3 states. This would ap- 
pear to suggest that one cannot rely upon discretisation 
corrections to continuum J = 1 operators to reliably pro- 
duce J = 3 states. In column (g) the continuum J = 4 
operator is also discarded and the J = 4 state vanishes. 
Finally, in column (h), the full operator basis, less the 
continuum J = 4 operator is used yielding a spectrum 
that lacks only the J = 4 state. 

In summary one should be sure to have operators with 
continuum overlap on to all the spins you expect to see. 
There is apparently little benefit, in terms of improving 
the precision of determination of the low-lying spectrum, 
of reducing the operator basis size by discarding opera- 
tors. 



C. Number of distillation vectors 

The results presented so far came from analysis of cor- 
relators computed on 16 3 lattices using 64 distillation 
vectors. We might wonder how the determination of the 
spectrum varies if one reduces the number of distillation 
vectors and thus reduces the computational cost of the 
calculation. This is particularly important given that, 
as shown in [5D], to get the same smearing operator on 
larger volumes one must scale up the number of distil- 
lation vectors by a factor equal to the ratio of spatial 
volumes. To scale up to a 32 3 lattice this would require 
64 x (||) = 512 vectors which is not currently a realiz- 
able number without using stochastic estimation |38j . 

In figure [9] we show the low- lying part of the extracted 
T l spectrum on the 743 16 3 lattice as a function of 
the number of distillation vectors used in the correlator 
construction. It is clear that the spectrum is reasonably 



stable for N > 32 but that the spectrum quality degrades 
rapidly for fewer vectors. 

It is also possible within distillation to implement a 
smearing operator other than □ = X}^=i£n£n- This 
is particularly relevant for large N where this smearing 
choice tends toward the identity and hence does not actu- 
ally filter out high-energy modes. An alternative smear- 
ing includes a weight function as n w — Y^,n=i w (^n)£n£b 
where w(X n ) might, for example, be a gaussian damping 

e~ CT A ™/ 4 . In figurej^we also show the spectrum obtained 
using N — 64 and a — 5.7 which, as shown in the inset, is 
a smearing radius that crudely approximates using only 
24 distillation vectors. The thus extracted spectrum dif- 
fers very little from the N = 64 spectrum, suggesting 
that with 64 vectors we are still far away from the "un- 
smeared" limit. 

One place (not shown in figure [9]) where the effect of 
reduction of the number of distillation vectors is seen 
clearly is for high-spin states. In particular the 4 
state seen at m/rriQ ~ 1.6 in T x is not reliably ex- 
tracted for N < 48. The need to have large num- 
bers of distillation vectors to reliably extract high spin 
mesons can be described in a simple free-field picture: 
In the continuum without gauge-fields, the eigenvec- 
tors of the laplacian (— V 2 £ = k 2 £) can be expressed 

as £(f) = = 4nE e i e 3i(kr)J2 m Y e m *(8,cj))Y e m (k)- 

From this expression it would appear that all £ values 
should contribute for any value of the eigenvalue k 2 and 
that high-spin states should be excited even by low eigen- 
vectors. But this argument does not take account of the 
radial behaviour which must be compared to the typical 
size of hadrons. The spherical Bessel function je(kr) is 
peaked at low values of kr with the peak position moving 
out to larger kr for larger I. Thus, since hadrons are only 
of finite size, i?, in order that the peak of ji(kr) remain 
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within r < R as £ increases, one must also increase the 
value of k. Thus, to have considerable amplitude multi- 
plying Y e and hence large overlap onto high-spin mesons, 
as £ increases one must include eigenvectors of higher 
eigenvalue, fc. 3 

In summary one is limited as to how few distillation 
vectors can be used if one requires reliable extraction 
of high-spin states. The results shown here suggest 48 
distillation vectors on a 16 3 lattice is the minimum, so 
96 distillation vectors on a 20 3 lattice are likely to be 
required. 



VIII. RESULTS 

As outlined in Section [TTJ our first results are obtained 
on lattices with pion masses between 400 and 700 MeV. 
The heaviest pion mass corresponds to the three-flavour 
symmetric point, while the lower masses have lower light- 
quark masses but the same bare strange quark mass; 
these are found to have only mild SU(3)f breaking, as 
indicated by mx/Tn^ which reaches only 1.4 on the light- 
est lattice, well below the physical value of 3.5. For each 
lattice (except one) we compute correlators on two vol- 
umes, 16 3 x 128 and 20 3 x 128. More complete details of 
the number of configurations, time sources and distilla- 
tion vectors used are given in Table [I] 



A. ~ 700 MeV results (743) 

This is an example of an exact S'?7(3)-flavour symmet- 
ric calculation. Since we are computing only the con- 
nected two-point correlators we will obtain the mass spec- 
trum of degenerate octets (e.g. degenerate pions, kaons 
and 77s). Results on 16 3 lattices have already been re- 
ported in Ref. 8J. To obtain singlet states like the 771 , 
which with exact flavour symmetry cannot mix with the 
rja, we would need disconnected two-point correlators. 

Variational analysis as described in previous sections 
leads to the irrep spectra from 16 3 lattices shown in fig- 
10 The colour-coding, tabulated in table|Vj indicates 
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the contin uum spin as determined by methods described 

overlaps are shown in figures [3] 



m section 



VI 



_ The A 

and|4j If the spin of a state is not unambiguously deter- 
mined by the methods of section VI it is represented by 



an orange box. On the other hand, if a state is extracted 
in the variational analysis but its mass cannot be accu- 
rately determined from fits to the principal correlator, it 
is represented by a grey box. 

In figure [TT| we show side- by-side the spectra obtained 
from 16 3 and 20 3 lattices. We note that there is really 
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TABLE V: Colour-coding used in spectrum irrep plots. 
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FIG. 10: Extracted spectra by irrep from 743 16 3 lattices. 



no change significantly outside statistical fluctuations be- 
tween the two volumes. We also note in passing that the 
128 distillation vectors used on the 20 3 lattices give a 
smearing that is essentially equivalent to that obtained 
by using 64 distillation vectors on the 16 3 lattices. In 
fact we observe that the independently extracted Z val- 
ues between the two volumes scale rather accurately as 



3 this analysis is somewhat oversimplified since in fact the size of 
hadrons is likely to increase with increasing spin 



Given the success of spin-identification we can sum- 
marise the results in a spectru m la belled by continuum 
J PC quantum numbers, figure 12 Here we show only 
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FIG. 11: Extracted spectra by irrep from 743 16 3 and 20 3 lattices. 



well-determined low-lying states. 

There are a number of notable features within this 
spectrum. Firstly there appears to be much of the 
n 2S+1 Lj distribution of non-exotic states predicted by 
the qq quark model (e.g. [39]). The left-hand pane of 
figure [12] has candidates for a ground-state S-wave pair 
(0 1 ) and a radial excitation at around m/rriQ ~ 
1.1. There is a complete (1,2,3) ,2 h D-wave set 
around m/ran ~ 1.2. The region around m/mo ~ 1.6 
appears to contain parts of a D-wave radial excita- 
tion and a G-wave ((3,4,5) ,4 h ) 4 . In the middle 
pane there is a clear P-wave, ((0, 1, 2) ++ , 1 H ) around 
m/mfi ~ 0.9 and a probable radial excitation near 
m/rriQ ~ 1.3. At m/niQ ~ 1.4 there is a complete -F- 
wave, (2, 3, 4)++, 3+-. 

In the right pane we see clearly a set of exotic J PC 
states, not accessible to a simple qq pair. Such states 
can be described in terms of constituents if additional 
degrees-of-freedom, either gluonic or extra qq pairs are 
included. If the extra component is gluonic the states 
are known as hybrid mesons and within models of glu- 
onic excitation, such as the flux-tube model [ID] , there is 
usually a roughly degenerate set of 1 h , H and 2" 1 
states. A hybrid nature for these states is suggested by 
the large overlap onto operators with essential non-trivial 
gluonic structure as described below. In this calculation 
we observe a 1 h lightest, with a pair of 0^ and 2" 1 
states nearly degenerate a little higher in mass. A second 
2+~ is then close to a second 1 h . States with exotic 
and 3 ^ quantum numbers are found to be considerably 
heavier, well above m/rriQ = 1.6. 



4 we have no operators capable of producing a spin-5 meson in the 
continuum. 



There are also a number of non-exotic quantum num- 
bered states which do not appear to fit into the qq 
n 2S+1 Lj classification. The ^,1 pair at m/rriQ ~ 
1.3 is probably too light to be the second radial excitation 
of the 5-wave (most likely the pair near m/mn ~ 1.55) 
and is partnered with a totally unexpected 2 h state. 
There may be excess states too in the positive par- 
ity sector but the situation is not totally clear above 
m/msj ~ 1.4. We note that the mass scale of these first 
excess states is comparable to that of the lowest lying 
exotic states shown in the right-hand pane. Further- 
more, these states have characteristically different over- 
lap behaviour compared to most other extracted states: 
they all overlap considerably onto operators featuring the 
commutator of two derivatives, that is the gluonic field- 
strength tensor. 

Given the large overlap onto operators requiring a non- 
trivial gluonic field distribution, we identify these states 
as hybrid mesons with non-exotic quantum numbers. 
Such non-exotic hybrids are predicted in models that as- 
sume non-trivial gluonic field configurations like the flux- 
tube model or constituent gluon models. In principle, 
such states can mix strongly with regular quark-model 
qq states leaving a spectrum which is not simple to in- 
terpret. In our results such large mixing may be present 
for the "excess" ^ state which has a large overlap also 
onto operators like ip^ij); however, this mixing does not 
appear to be present to the same degree for the 1 state. 

Detailed model-dependent interpretation of the spec- 
trum, comparing the overlap values with the expectations 
of a bound-state quark model (analogous to that done for 
charmonium in [41] ) and considering the degree of mixing 
of non-exotic hybrids and quark model states will follow 
in a subsequent publication. 
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FIG. 12: Spin-identified spectrum of isovector (octet) mesons from the 743 lattices. 16 3 (solid) and 20 3 (dashed) spectra agree 
well. Ellipses indicate that there are heavier states with a given J PC but that they are not well determined in this calculation. 



B. Quark mass dependence 

Here we move away from the SU(3) flavour point by 
lowering the mass of two degenerate "light" flavours and 
keeping one remaining strange flavour heavy. We have 
access to isovector mesons from the connected correla- 
tors with a light quark and a light antiquark, and kaons 
from the connected correlators with a light quark and a 
strange antiquark. It is also possible for us to compute 
the connected part of correlators with both quark and an- 
tiquark being strange, the so-called "strangeonium" . We 
recognise that neglecting the disconnected contributions 
to these diagrams leads to a non-unitary description of 
this particular element of our calculation. Of course it 
is also true that "strangeonium" states are not necessar- 
ily QCD eigenstates since being isoscalars they can mix 
with light-light isoscalars and pure-glue states through 
disconnected diagrams 5 . The classic extreme examples 
are the rj, v( system which is mixed almost as SU (3) p 
octet-singlet and the w, </> system which is mixed almost 
as 1 1, ss. 

In the figures we show extracted state masses as a 
function of £n = | |°'™^j 2 which we use as a proxy for 
the quark mass [29 . The state masses are presented via 



5 work is ongoing within the Hadron Spectrum Collaboration |20| 
to utilise distillation methods to efficiently compute disconnected 
two-point functions, allowing extraction of the true QCD eigen- 
states 



«m7 s \ The ratio of the state mass (m H ) to the fi- 
baryon mass computed on the same lattice removes the 
explicit scale dependence and multiplying by the phys- 
ical 57-baryon mass conveniently expresses the result in 
MeV units. This is clearly not a unique scale-setting 
prescription, but it serves to display the data in a rela- 
tively straightforward way. We remind the reader that 
the data between different volumes and quark masses are 
uncorrelated since they follow from computations on in- 
dependently generated dynamical gauge-fields. 



1. Isovector mesons 

Figures fl3| [14} [T5| [T6| show the extracted spin- assigned 
spectra for mesons of isospin-1 having a range of J PC 
quantum numbers (the neutral members of the 1=1 
triplet are eigenstates of C, the charged members are 
eigenstates of G-parity with G = — C). We use the PDG 
[4*2"] nomenclature for meson states throughout. In those 
cases where two states are almost degenerate we shift one 
in the horizontal direction by an amount S£q = 0.005 for 
clarity. In some cases, for comparison, we plot the mass 
of the lightest meson-meson pair which in S'-wave would 
have the appropriate quantum numbers - the mass follows 
from the simple sum of the extracted masses on these 
lattices. This may involve a so-far undetermined isoscalar 
mass and in these cases we use the approximations m u rj 
m p , rrif 2 w m a2 and the crudely estimated r\ mass from 
[29] , Occasionally we extract a low- lying state that is 
reasonably robust against the changes in analysis detailed 
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in section VII but whose principal correlator is quite 
noisy leading to a relatively poorly determined mass - we 
show these states with dashed symbols. 

Not shown are results for isovectors, the poi which 
are exotic. The lightest such state we extract is at least 
2 GeV heavier than the p at all our quark masses. The 
exotic 3 % the 7r 3 , is found to be similarly heavy. 



2. Kaons 



In the kaon sector, we no longer have charge- 
conjugation as a good quantum number, with only J p 
remaining in the continuum which is then subduced into 
on a cubic lattice. We compute a correlator matrix 
for a given A p using the concatenated list of all A p+ and 
A p ~ operators. 

Using a combination of experiment and models 
there are suggestions that resonant kaon states are 
mixtures of basis states of opposite C with a rather 
large mixing angle. For example the axial kaons, 
i£i (1270), i£i (1400), are suggested to be mixtures of ba- 
sis states K\a(C = +), K\b(C = — ) with a mixing angle 
near 45°. Clearly this mixing relies upon being far away 
(rn^ ys /mP hys = 3.5) from the SU(3) p limit, since in that 
limit there is effectively restoration of (a generalisation 
of) C-symmetry. All the lattices presented in this paper 
can be considered to be rather close to the SU(3)p limit 
(1 < TOjf/m w < 1.39) and we observe little or no mixing. 
This is suggested at the correlator level (see figure 17 for 
the 840 16 3 correlator), and verified in the Z values of 
the spectrum extraction (see figure 18 for the 808 16 3 
lattice (left) and the 840 16 3 lattice 



for the 
right)). While the 



840 lattice shows a greater degree of opposite C mixing 
than the 808, indicating an increased breaking of SU(3) f 
symmetry, the mixing is still very small in absolute terms 
and states are approximately eigenstates of C. The 840 



kaon spectrum is shown in figure 19 with the dominant 
C-eigenstate noted for each state. 

The light-quark mass dependence of kaon states is dis- 
played in figures [20 21 22 where color-coding indicates 
the dominance of C-eigenstates within the spectrum. 



In figure [23] we show only a subset of possible 
"strangeonium" J PC , those for which there is some phe- 
nomenological evidence that the QCD eigenstates are in 
fact close to being pure ss. We note that using this par- 
ticular scale-setting scheme we observe light-quark mass- 
dependencies that are very flat. 
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FIG. 17: Normalised corrrelation matrix (Ctj / \ZCuCjj) on 
timeslice 5 in the kaon irrep on 840 16 lattice. Operators 
with C = — first, C — + second. 
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FIG. 18: Overlaps, Z, of kaon operators onto lightest two 
J = 1, lightest two J — 3 and lightest J = 4 states for 808 
16 3 (left) and 840 16 3 (right) lattices. Operators with C = ± 
are grouped to show clear separation. Normalisation as in 
figure [3] 
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FIG. 20: Lightest negative parity kaons. Color coding indicates dominance of a particular charge-conjugation eigenstate. 
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IX. MULTI-MESON STATES 

In the previous section we presented the extracted 
spectra from calculations with four different light quark 
masses on two different lattice volumes. In each case we 
were able, using the operator overlaps, to match states 
across irreps that we believe are subduced from the same 
continuum spin state. This suggests an interpretation 
of the spectrum in terms of single-hadron states, while 
in principle our correlators should receive contributions 
from all eigenstates of finite- volume QCD having the ap- 
propriate quantum numbers. This includes multi-meson 
states which in finite volume have a discrete spectrum. In 
fact, in a theory of non-interacting mesons the spectrum 



m[A(p)B(-p)] = ^m\ + W + ^ml + \p\\ (5) 

where only discrete values of the momentum p are allowed 
by the boundary conditions: a s p = j^- (n x , n y , n z ^j (L s is 
the spatial lattice extent in lattice units, i.e. 16 or 20 for 
the lattices we are using). Clearly this spectrum, with 
the exception of the states with p = (0, 0, 0), will change 
considerably under changes in volume. Ongoing work 
by the Hadron Spectrum Collaboration is presented in 
Ref. [35] . 

Within QCD mesons interact and this interaction has 
a range of possible forms, for example: the repulsive in- 
teraction of two pions in the isospin 2 channel, the some- 
what attractive interaction of two pions in isospin that 
gives rise to the a enhancement and the strong attraction 
in isospin 1 that corresponds to the relatively narrow p 
resonance. As shown by Liischer[27j, taking account of 
hadron interactions, the finite volume energy spectrum 
will be modified with respect to equation [5] The mod- 
ification (at least in the elastic case) can be related to 
the hadron-hadron scattering phase shift which encodes 
details of attractive or repulsive interactions and even of 
resonant behaviour. 

A simple schematic framework in which to view the 
finite- volume eigenstates is in terms of admixtures of ide- 
alised non-interacting basis states. For example, at low 
energy in the T 1 channel, one might consider there 
to be a space of non-interacting pion-pair states with 
the various relative momenta allowed in a finite cubic 
box. In addition we can allow a space of single-hadron 
vector bound states like the p, which we assume to be 
localised to a region of space somewhat smaller than 
the size of the lattice box. The pion-pair state ener- 
gies vary rapidly with changing box-size while the "p" 
bound-states would be essentially volume-independent 
for volumes larger than the size of the bound-states. If 
one supplies a resonant phase shift (such as the p in 
7T7T scattering), Liischer's formulae give rise to avoided 
level crossings as a function of lattice size that resemble 
the behaviour of approximate eigenstates within time- 
independent quantum mechanical perturbation theory. 
The finite-volume eigenstates can be viewed then as ad- 
mixtures of the "p" bound-states and the pion-pair states 
where the degree of mixing is determined by the phase- 
shift and the volume of the box. 

In figures [24] and [25| we show the extracted 743 spec- 
trum superimposed with the positions of non-interacting 
meson-meson states for 16 3 and 20 3 lattices. The distri- 
bution across irreps is determined using the "in-flight" 
cubic symmetry group theory tables from [JS] |3S] ■ We 
show only pairs of SU(3)f octet states (since we have 



6 we are neglecting small discretisation effects in the dispersion 
relation. 
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not determined the masses of the singlets or any other 
possible multiplets) and do not indicate the multiplicity 
of flavoured states for each level (which follows from the 
SU(3) Clebsch-Gordan series for 8 ® 8 -> 8). 



What is clear from Figures 24 and 25 is that the ex- 
tracted spectrum does not seem to be related in any ob- 
vious way to the non-interacting two-particle spectrum. 
The distribution of two-particle states across different ir- 
reps is determined not by the cubic symmetry of our 
discretized lattice, but rather by the spatial momenta al- 
lowed by the boundary conditions on the cubic volume 
in which we are performing our calculations. We have 
seen that the observed spectrum split across the differ- 
ent irreducible representations conforms to that expected 
for single-particle states, with only negligible effects from 
cubic symmetry on the scale a s . In contrast, the pattern 
across different irreps expected for multi-particle states 
would be quite different. This leads further credence to 
our assertion that two-particle states are contributing lit- 
tle to the calculated correlators. 

How can we explain the lack of two-meson states in 
the spectrum we observe? Clearly, the basis set of op- 
erators employed has a very small overlap with states 
that predominantly resemble a composite of two mesons 
with well-defined and opposite momentum. This small 
overlap means the contribution to the correlation func- 
tion coming from the significant number of two-meson 
states is not resolved within the statistical precision of 
our calculation. 

This effect has been observed before in other 
dynamical-quark Monte Carlo measurements [47l [48] of 
the energies of states that are above threshold, although 
perhaps not to the dramatic degree seen in this study. 
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FIG. 25: As figure 24 for A+ 



One cause of this supression may come directly from con- 
finement dynamics of strongly interacting gluons. The 
inter-quark gluon flux shows considerable reluctance to 
break in the regime where it is energetically favourable 
to form two well-separated color singlets. This has been 
seen |49l I50j as a very small overlap of an operator com- 
prising a gluonic string onto the ground-state of this sys- 
tem. 

The overlap of a localised quark-bilinear operator onto 
a two-meson state will be supressed by l/y/V, where V 
is the lattice volume, if the operator creates a resonance 
with a finite width in the infinite volume limit. This fall- 
off is matched by a growth in the density of states with 
the volume and the resonant state thus maintains a finite 
width as the mixing with each discrete state falls. The 
simulations in this study are carried out in cubic vol- 
umes with side-lengths bigger than 2 fm, which might be 
sufficiently large that the mixing between one of the low- 
lying two-particle states and a resonance is suppressed 
sufficiently for it to be undetectable with the quark bi- 
linear operator basis. 

Even if the mixing between localised single-hadron 
states and two-meson states to form resonance-like finite 
volume eigenstates is not small, there still remains a prac- 
tical difficulty associated with using only quark bilinear 
operators. In this case the state can be produced at the 
source timeslice through its localised single-hadron com- 
ponent, while the correlator time dependence obtained 
from e~ Ht will indicate the mass of the resonant eigen- 
state. Consider a hypothetical situation in which a sin- 
gle two-meson state, denoted by |2), mixes arbitrarily 
strongly with a single localised single-hadron state, |1), 
with all other states being sufficiently distant in energy 
as to be negligible. There will be two eigenstates 

| a) = cos0|l) +sin#|2) 
\b) = -sin0|l) + cos0|2>, 
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with masses m a ,mt,. At the source (and sink) only the 
localised single-hadron component of each state overlaps 
with the operators in our basis and hence the overlaps, 
Z?' = (a, b|Oi|0), will differ only by an overall multi- 
plicative constant, Z? = cos0Zf\ Z\ = -sm9Z\ 1} . As 
such the eigenvectors v a ,v point in the same direction 
and cannot be made orthogonal. Thus the time depen- 
dence of both states will appear in the same principal 
correlator as 

X(t) ~ A a e~ m ' {t ~ to) + A b e- mb{t - to) + ... 

Since m a and m& most likely do not differ significantly 
(on the scale of a^ 1 ) it will prove very difficult to ex- 
tract a clear signal of two-exponential behavior from the 
principal correlator. This is precisely why the varia- 
tional method's orthogonality condition on near degen- 
erate states is so useful, but we see that it cannot work 
here and we are left trying to extract two nearby states 
from a x 2 fit to time-dependence. Typically this is not 
possible and reasonable looking fits to data are obtained 
with just one low-mass exponential. 

If this is what is occurring in parts of our extracted 
spectrum, that we are extracting states which are ad- 
mixtures of "single-particle" and multi-meson states, but 
that we are not extracting the orthogonal mixtures, how 
should we interpret the mass values we are extracting? 
One conservative approach would be to suggest that our 
mass values are accurate only up to the hadronic width 
of the state extracted, since this width is correlated with 
mixing with multi-meson states via the scattering phase- 
shift. 

We would like to explicitly observe resonant behaviour 
of states in our calculations and as such we must coun- 
tenance the inclusion in our operator basis of operators 
with larger numbers of fermion fields in order to obtain 
healthy overlap with multi-hadron states. This can be 
done while respecting the lattice symmetries using the 
tables in [45J 06] . By using single-meson operators sub- 
duced into the 'in flight' little-group irreps from opera- 
tors of definite continuum helicity, it may prove possi- 
ble to utilise something similar to the spin-identification 
carried out in this paper. These constructions are under- 
way, and distillation, with the possible use of a stochastic 
estimator [35], affords an efficient numerical means of im- 
plementing them [32]. 



X. SUMMARY 

We have described in some detail our method for ex- 
tracting a large number of excited states using dynamical 
anisotropic lattices, distillation technology and a varia- 
tional analysis with an extensive basis of carefully con- 
structed operators. We have demonstrated the stability 
of the spectra with respect to changing the specific set of 
meson operators used, the number of distillation vectors 
and the details of the variational analysis. Our method 



of spin identification based on operator overlaps has en- 
abled us to confidently assign continuum spin to many 
states. 

We have successfully applied these techniques for two 
volumes on multiple mass sets: one with three degen- 
erate flavours of quark (743) and three with lighter u 
and d quarks giving mild breaking of SU(3)f and pion 
masses down to ~ 400 MeV. In all cases we see mostly 
no significant volume dependence. On all mass sets and 
volumes we are able to reliably extract a large number of 
excited states with all PC combinations, states with high 
spin, up to and including spin four (4 ++ , 4~ + and 4 ), 
and states with exotic quantum numbers (0 -1 , I *~ and 
2 H ). The exotic states are particularly interesting and 
their presence points to the influence of explicit gluonic 
degrees of freedom. In Fig. [26] we summarise our results 
on exotic states and compare with previous lattice QCD 
results from Refs. [2IH2S1- 

The extracted spectra show features of the n 2S+1 Lj 
state assignment of bound-state quark models, along 
with states (both exotic and non-exotic) which do not 
seem to lie within that classification. A detailed model- 
dependent interpretation of these spectra is called for, 
comparing overlaps with quark model expectations and 
determining the degree of mixing of non-exotic hybrids 
and quark model states. This work is ongoing. 

We have presented kaon spectra and observe little or no 
mixing between the two charge conjugation eigenstates 
(C = + and C = — ); the resulting spectrum largely 
corresponds to the superposition of the C = — and C = 
+ isovector spectra modulo the mass shift due to the light 
quark - strange quark mass difference. In the SU(S) p 
limit there can be no such mixing and on the mass sets 
considered we are still rather close to this limit (f < 
mx/mir < 1.39), so it is therefore not surprising that the 
mixing is small. Of particular interest at lighter quark 
masses, closer to the physical m^/m^ — 3.5, will be the 
mixing between the axial kaons, l£"i(1270) and ifi(1400), 
which, using a combination of experiment and models [4"3l 
|4*4"] is expected to be large. 

We have argued that we see little evidence for two- 
particle states in our spectra and that to study such 
states we need to construct operators with a larger num- 
ber of fermion fields. Such constructions are in progress 
and we believe that the addition of these operators 
will lead to a denser spectrum of states which can be 
interpreted in terms of resonances via techniques like 
Liischer's and its inelastic extensions. 

With the excited state spectra extracted herein, we ar- 
gue that it does not make sense to attempt chiral extrap- 
olation given that we cannot form a clear field-theoretic 
interpretion of the extracted energy levels. Once we have 
a handle on the two-meson levels, we can apply the tech- 
niques mentioned above to extract something like a phase 
shift, or more generally elements of the S'-matrix, at dis- 
crete energy values. The phase shift may show resonant 
behaviour, which can be fitted with, in the simplest case, 
a Breit-Wigner form. The mass and width parameters 
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of this Breit-Wigner are quantities which should be more 
amenable to chiral extrapolation. 

A further avenue of study is the computation of dis- 
connected two-point correlators giving access to isoscalar 
mesons; here we are interested in determining how QCD 
decides to mix light and strange. In addition, methods 
similar to those detailed in this paper are being applied in 
the baryon sector where the lattice irrep spectrum suffers 
from an even greater degree of degeneracy. An important 
aim of the Hadron Spectrum Collaboration is the calcu- 
lation of light meson photocouplings which are relevant 
for, amongst other things, the GlueX experiment at the 
JLab 12 GeV upgrade where light mesons will be studied 
in photoproduction, with particular interest in exotics. 
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Appendix A: Subduction Coefficients 

Here we give a derivation of the subduction coefficients 
before listing their explicit values for all integer spins up 
to J = 4. An alternative derivation using the group 
theoretic projection formula is also described. 

The continuum spin J is reducible under the group 
of lattice rotations (the octahedral group or equivalcntly 
the cubic group). We use "subduction" coefficients to 
project the continuum based operators onto their suitable 
octahedral group based versions via 



O 



[J] 

A, A 



(Al) 



M 



where J ' M are the continuum operators with some defi- 
nite total spin J and spin component M . For each J — > A 
there is a matrix in the values of M and the rows of the 
irrep, A, that performs this mapping, i.e. the subduction 



coefficients, S 



J,M 
A, A • 



As shown in Table III the J = continuum spin sub- 



duces only onto the one-dimensional Ai irrep and so triv- 



ially we have S 



0,0 



1. The simplest non-trivial subduc- 



tion is that of continuum J ■■ 
is given by 



(•/=!)-» T a 



1 into the T\ irrep and this 
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(A2) 



where the T\ is a faithful representation of J = 1 (The 
basis used here follows Ref. [2])- 



The subduction coefficients for higher spins can be built up using continuum and octahedral group Clebsch-Gordan 
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coefficients using 



S 



J,M 
A, A 



Ai,A 2 "»i,"i2 \ 



A Ai A 2 
A 2 , 



or Jl qJ 2 ,m2 ; 



'Ai,Ai -A 2 ,A 2 ! ( J i' m i; J 2 ,m 2 |J,M) . 
Here ( Ji, Mi; J 2 , M 2 | J, M) is the usual SO (3) Clebsch-Gordan coefficient for Ji <£> J 2 J and C 



(A3) 



is the 



A Ai A 2 
A Ai A 2 

octahedral group Clebsch-Gordan coefficient for Ai <g> A 2 — » A. iV is a normalisation factor, fixed by the requirement 

5a,a'- 

-.J.M 



that the subduction coefficients form an orthogonal matrix, Em^a'a ^A' a* = ^A,A'<^. 



M A, A "A', A' 

This iteration formula can be constructed by noting that, for appropriately normalised states, S^a = M^> M) 
and so 



S A,x= N Yl E (A,A|Ai,Ai;A 2 ,A 2 )(Ai,Ai;A 2 ,A 2 |Ji,mi;J 2 ,m 2 )(Ji,TO i; J 2 ,m 2 |J,M). (A4) 

Ai ,A 2 mi ,m2 

Substituting for the continuum and octahedral group Clebsch-Gordan coeffi cients and (Ai, Ai; A 2 , A 2 | J\, mi; J 2 , m 2 ) = 
(Ai, Ai| Ji,mi) (A 2 , A 2 | J 2 , m 2 ) = S^'^S^'™* gives the result in Eq. 
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For J = 2 to J = 4 the subduction coefficients are shown below: 
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An alternative method for constructing the subduction coefficients is by using the group theoretic projection formula: 

,J,M' 



QG — : ^ — ' 



(A5) 



2G 



where G is the octahdedral group, go is the order of the group G (i.e. 24), c?a is the dimension of octahedral group 
irrep A, Rmm' is an element of G acting on J ' M and Tj^ (R) is the representation of R in A. Here the operators 

O ' have definite spin and so Rmm> = D$, M (R) (a Wigner-£> matrix). The T^(R) are found, for example, by 
considering a basis for the irreps in terms of spherical harmonics [5] and then using the transformation properties of 

spherical harmonics under rotations. Once all possible operators have been found (considering all M and A), the 
linearly independent combinations are constructed. These combinations then give the subduction coefficients which 
are identical, up to possible phases and choice of basis, to t hose obtained using the method described above. 



